/*
 *  lib-smv: Matrix - Vector Library targeted to Structural Mechanics Problems
 *  Copyright (C) 2006-2008  St.Brcic, Lj.Zugic
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
/** @file Cholesky.h
 *  @brief Definition of class Cholesky - LU decomposition of SPD matrices
 *
 *        Standard decomposition     A = L * LT 
 *        Generalized decomposition  A = L * D * LT
 *
 *  DIRECT SOLUTION OF LINEAR EQUATIONS A * x = b or A * X = B
 *  
 *  LU decomposition of symmetric positive definite (SPD) matrices used for:
 *  (1) Solution of linear equations
 *  (2) Inverse of a SPD matrix
 *  (3) Determinant of a SPD matrix
 * 
 */ 

#ifndef _CHOLESKY_H
#define _CHOLESKY_H

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <cassert>
#include <string>
#include "Arrays.h"
#include "iArrays.h"

using namespace std;

namespace smv_lib
{

//! standard or generalized Cholesky decomposition
enum DECOMPOSE {chol_STANDARD, chol_GENERAL};
//------------------------------------------------------------------
/**
** \brief class Cholesky - LU decomposition of symmetric matrices
**        Standard decomposition     A = L * LT 
**        Generalized decomposition  A = L * D * LT
**
**  DIRECT SOLUTION OF LINEAR EQUATIONS A * x = b or A * X = B
**  
**  LU decomposition of symmetric positive definite (SPD) matrices used for:
**  (1) Solution of linear equations
**  (2) Inverse of a SPD matrix
**  (3) Determinant of a SPD matrix
**  
**/
//-------------------------------------------------------------------
class Cholesky
{
  private:
    //! Matrix A (should be square real and positive definite)
    Matrix A;
    //! Order of matrix A
    int N;
    //! Vector b - right-hand side of equations A * x = b
    Vector b;
    //! Vector x - solution of linear equations A * x = b
    Vector x;
    //! Number of right-hand side vectors (no. of eqs. with the same coefficient matrix A)
    int M;
    //! Matrix B - right hand side of equations A * X = B (B is matrix N x M)
    Matrix B;
    //! Matrix X - the solution of linear equations A * X = B (X is matrix N x M)
    Matrix X;
    //! Inverse of the matrix A
    Matrix* A_inv;
    //! Determinant of the matrix A
    double Det;
    //! Lower triangular Cholesky factor (A=L1*L1T or A=L2*D*L2T)
    Matrix L1, L2;
    //! Transpose of lower triangular Cholesky factor (A=L1*L1T or A=L2*D*L2T)
    Matrix* LT1, *LT2;
    //! Vector with diagonal elements in A = L * D * LT
    Vector D1, D2;
    //! Inverse of lower triangular factor
    Matrix* L_inv;
    //! Indicator: true if matrix A is given
    bool isMatrix_A;
    //! indicator: true if matrix A is symmetric
    bool isSymmetric;
    //! indicator: true if matrix A is positive definite
    bool isPosDefinite;
    //! indicator: true if matrix A is SPD (symmetric and positive definite)
    bool is_SPD;
    //! Indicator: true if vector b is given (A*x=b)
    bool isVector_b;
    //! Indicator: true if matrix B is given (A*X=B)
    bool isMatrix_B;
    //! Indicator: true if more sistems of eqs. with the same matrix A are solved
    bool isMore;
    //! output stream
    ofstream out;
    //! Name of the problem (output file is Name + ".txt")
    string Name;
    //! Indicator: true if the name of the problem is given
    bool isNameGiven;
    //! Standard (A=L*LT) or generalized decomposition (A=L*D*LT)
    DECOMPOSE Decompose;
    //! Indicator: true if standard decomposition: A = L * LT
    bool isStandard;
    //! Indicator: true if generalized decomposition: A = L * D * LT
    bool isGeneralized;
    //! Indicator: true if decomposition is over
    bool isDecomposed;
    //! Indicator: true if inverse of the Matrix A is obtained
    bool isInverse;
    //! Indicator: true if linear equations are solved
    bool isSolved;

  public:
    //! empty (default) constructor
    Cholesky(void);
    //! constructor with a given matrix A and decomposition method (standard is default)
    Cholesky(const Matrix& AA, DECOMPOSE dec = chol_STANDARD);
    //! constructor with a given matrix A and vector b - to solve equation A*x=b
    Cholesky(const Matrix& AA, const Vector& bb);
    //! constructor with a given matrix A and matrix B - to solve equations A*X=B
    Cholesky(const Matrix& AA, const Matrix& BB);
    //! destructor
    ~Cholesky(void);

    //----------------------------- set and get functions
    //! set Matrix A and decomposition method (standard is default)
    void setMatrix_A(const Matrix& AA, DECOMPOSE dec = chol_STANDARD);
    //! set the method of decomposition: standard or generalized
    void setMethod(DECOMPOSE dec);
    //! set the name
    void setName(const string name) {
      Name = name;
      isNameGiven = true;
    }
    //! set the right-hand side vector b
    void setVector_b(const Vector& bb);
    //! set the right-hand side matrix B
    void setMatrix_B(const Matrix& BB);
    //! get the method of decomposition: standard or generalized
    DECOMPOSE getMethod() {
      return Decompose;
    }
    //! get lower triangular factor L (either A=L*LT or A=L*D*LT)
    Matrix get_L();
    //! get upper triangular factor L (either A=L*LT or A=L*D*LT)
    Matrix get_LT();
    //! get diagonal elements of decomposition (for A=L*D*LT)
    Vector get_D();
    //! get inverse of lower factor (either A=L*LT or A=L*D*LT)
    Matrix get_L_inv();
    //! get inverse of the matrix A
    Matrix getInverse();
    //! get determinant of the matrix A
    const double getDeterminant() const;
    //! get solution of linear equations A * x = b
    Vector getSolution();
    //! get one solutions of the system of linear equations A * X = B
    Vector getSolutionNo(int i);
    //! get all solutions of the system of linear equations A * X = B
    Matrix getAllSolutions();
    //! give the method of decomposition
    const string giveMethod() const;

    //---------------------------------- working functions
    //! decomposition of the matrix A
    void Decomposition();
    //! inverse of the matrix A
    void Inverse();
    //! solution of linear equations A * x = b or A * X = B
    void Solution();
    //! write results to predefined file "CHOLESKY.txt"
    void Write();
    //! write results to file with a given name
    void WriteTo(const char* file);

  private:
    //! standard Cholesky decomposition
    void cholDec();
    //! generalized Cholesky decomposition
    void cholDecGen();
    //! inverse matrix
    Matrix cholInv();
    //! solution of equations A * x = b
    Vector cholEqs(const Vector& bb);
    //! solution of system of linear equations A * X = B
    Matrix cholEqs(const Matrix& BB);
    //! determinant of the matrix A
    const double Determinant();
    //! true if the input is OK
    bool checkInput();
    //! true if matrix is symmetric
    bool checkSymmetric();
    //! true if symmetric matrix is positive definite
    bool checkPosDefinite();
    //! control of Cholesky decomposition
    Matrix cholControl();
};

}   // end of namespace smv_lib

#endif  // _CHOLESKY_H

//============================== end of file

